TUTORIAL 7 : Introduction to lattice-switch Monte Carlo

Author:Tom L. Underwood, t.l.underwood{at}bath.ac.uk

Introduction

A common problem is as follows: given two solid phases, which is more stable at a given temperature and pressure? The more stable phase is that with the lower free energy of the two. Thus calculating the free energy difference between the two phases reveals which phase is more stable. Lattice-switch Monte Carlo (LSMC) is a powerful method for calculating such free energy differences. Here we demonstrate how to use LSMC in DL_MONTE to calculate the free energy difference between two solid phases at a given temperature and pressure.

Specifically, we use LSMC to investigate which of the fcc or hcp crystal phases of the Lennard-Jones system are stable on the \(P=0\) isobar. The question of whether the fcc or hcp phase is the stable phase of the Lennard-Jones solid has a long history, and was resolved only relatively recently. The answer is that the hcp phase is stable at lower temperatures and pressures, while the fcc phase is stable at higher temperatures and pressures (below the melting line). We will use LSMC to investigate whether this behaviour is borne out in simulations using 216-particle cell and the following incarnation of the Lennard-Jones system: each particle in the cell interacts with all other particles via the Lennard-Jones potential

\[\phi(r) = 4\epsilon\biggl[\Bigl(\frac{\sigma}{r}\Bigr)^{12}-\Bigl(\frac{\sigma}{r}\Bigr)^{6}\biggr]\]

with no cut-off, using the minimum image convention (and there is no other contribution to the total energy, e.g. due to tail corrections).

Prerequisites

This tutorial assumes that the reader is familiar with the basics of DL_MONTE, i.e. can perform simple NVT and NPT simulations. Moreover it is also assumed that the reader has a basic understanding of biased sampling techniques.

Solutions

Note that example output and input files for all simulations to be performed in this tutorial can be found in the solutions directory. This is important in light of the fact that some of the simulations could take a long time to run: one may wish to skip the actual running of the simulation and proceed straight to the analysis or inspection of the output files.

Background - LSMC

We begin by providing a short description of LSMC. More detailed information can be found in the user manual, as well as the references provided at the end of the tutorial.

Consider the free energy difference between two phases 1 and 2. This is given by the following equation:

\[\Delta G \equiv G_1 - G_2 = -k_BT\log(p_1/p_2)\]

where \(p_1\) is the probability being in phase 1 and \(p_2\) is the probability of the system being in phase 2 (assuming that the system can only be in either phase 1 or phase 2); \(G_1\) is the free energy of phase 1, and similarly for \(G_2\); and \(T\) is the temperature of the system. Note that we are treating a system at constant temperature and constant pressure here, so the free energy here is the Gibbs free energy. In principle the above equation could be exploited to calculate \(\Delta G\) from simulation as follows: measure the times \(t_1\) and \(t_2\) that the system spends in phases 1 and 2 during the simulation, and then use the fact that \(t_1/t_2=p_1/p_2\) in the above equation to obtain \(\Delta G\). Of course this approach requires that the both phases are explored on timescales accesible to simulation. Unfortunately this is not the case in conventional methods which invoke ‘realistic’ particle dynamics. In such methods the time taken for the system to transition from one phase to the other and back again is simply too long for the above approach to be tractable.

The key aspect of LSMC is a lattice switch Monte Carlo move which, if accepted, takes the system directly from the current phase to the ‘other’ phase. In theory, if lattice switch moves are attempted and accepted frequently, then the system would explore both phases on reasonable timescale. This in turn would enable \(\Delta G\) to be evaluated via the above equation. In a lattice switch move a configuration corresponding to the other phase is constructed from the current configuration. The new configuration is constructed by ‘switching’ the underlying lattice which characterises the current phase for a lattice which characterises the other phase, while preserving the displacements of all particles from their lattice sites. This is best illustrated pictorially; see the figure below.

Illustration of a lattice switch move. The current configuration (left) corresponds to the square phase, and the move will take the system to a new configuration in the triangular phase (right). Note that in the current configuration the particles are all close to their lattice sites (red crosses, which form a square lattice). The move transforms the underlying square lattice (red crosses) into a triangular lattice (blue crosees), while keeping the displacements of the particles from their lattice sites (black arrows) unchanged - the displacement of particle \(n\) is the same in the current and new configurations.

However this is not the whole story. While lattice switch moves provide a means for generating reasonable-looking configurations in the other phase, these configurations are usually of very high energy, and accordingly the Metropolis algorithm almost always rejects lattice switch moves. Thus lattice switch moves alone are not enough to enable both phases to be explored in a single simulation. The extra ingredient which is required is biased sampling, which we exploit to increase the proportion of lattice switch moves which are accepted to an acceptable level. To elaborate, we first define an order parameter \(M\) which distinguishes the following three types of configuration:

1. configurations realised in phase 1 at equilibrium; those which one would see in a conventional Monte Carlo simulation of phase 1. These configurations have large negative values of \(M\), i.e. \(M\ll 0\). We refer to these configurations as phase-1 equilibrium configurations
2. configurations realised in phase 2 at equilibrium; those which one would see in a conventional Monte Carlo simulation of phase 2. These configurations have large positive values of \(M\), i.e. \(M\gg 0\). We refer to these configurations as phase-2 equilibrium configurations
3. configurations where lattice switch moves both from, and to the configuration are highly likely. As described in a moment, these configurations act as a gateway between the two phases, for which reason we refer to them as gateway configurations. They have small values of \(M\), i.e. \(M\approx 0\)

By defining \(M\) in this way we have created a pathway linking the two phases; traversing \(M\) takes us between phases. Consider a phase-1 equilibrium configuration. Here, \(M\ll 0\), and lattice switch moves are fruitless - they are essentially always rejected. Gradually increasing \(M\) to \(\approx 0\) however takes the system to gateway configurations where lattice switch moves are likely to be accepted, and the system can transition (via a lattice switch move) to phase 2. Increasing \(M\) further eventually takes us to the phase-2 equilibrium configurations, where lattice switch moves are again fruitless. With this in mind we perform a simulation using biasing to explore the whole range of \(M\) uniformly. The result is that both the equilibrium configurations for each phase, and the gateway configurations which are necessary for sucessful switching to and from both phases, are sampled. There is of course the question of how \(M\) is defined so that it has the desirable properties just mentioned. This is beyond the scope of the current tutorial; it is not esssential to know this to effectively use LSMC in DL_MONTE. However the references provided at the end of the tutorial for further information regarding how \(M\) is defined.

LSMC workflow

We now turn to the practicalities of LSMC free energy calculations. Calculating a free energy difference using LSMC in fact involves a number of steps:

1. A preliminary simulation for each phase. Above we said that we would use biased sampling to sample the whole range of \(M\) uniformly. Of course, in practice we must choose upper and lower limits to the range of \(M\) which our simulation will explore. However, we do not know a priori what limits are appropriate. We must therefore perform a short preliminary simulation for each phase to determine the appropriate range of \(M\) to use in subsequent simulations. These simulations should be ‘conventional’ Monte Carlo simulations, i.e. lattice switch moves and biasing should not be invoked. However the LSMC order parameter \(M\) must be tracked during the simulation.
2. A bias-generation simulation. Using biased sampling to sample uniformly over the range of \(M\) (determined in the preliminary simulations) requires knowledge of the bias function which results yields uniform sampling. Unfortunately this too is not known a priori. Accordingly the next step is to perform a bias-generation simulation to learn the bias function which leads to uniform sampling over our chosen range of \(M\).
3. A production simulation: a biased simulation using the bias function determined in the last step, which samples uniformly over \(M\) in the chosen range, and accordingly explores both phases. From the data obtained from this simulation we can determine \(p_1/p_2\), and hence \(\Delta G\) using the equation given at the beginning of this section.

Below we provide a walk-through of this procedure for calculating the hcp-fcc free energy difference of the Lennard-Jones system described at the beginning of this tutorial at temperature \(T=0.1\epsilon/k_B\) (where \(\epsilon\) is the Lennard-Jones well depth and \(k_B\) is Boltzmann’s constant) and pressure \(P=0\), using DL_MONTE.

Preliminary simulations

A basic LSMC simulation in DL_MONTE requires 5 input files: CONTROL, FIELD, CONFIG, CONFIG.1 and CONFIG.2. The first 3 of these files are required in all DL_MONTE simulations. However CONFIG.1 and CONFIG.2 are specific to LSMC. What follows is a description of the input files for our first preliminary simulation. These files can be found in the main tutorial directory.

FIELD

We begin with the FIELD file (note that the line numbers at the start of each line are not actually present in the file):

Lennard-Jones potential with sigma=epsilon=1, cut-off infinite, and 2 configs for PSMC, energy units in k_BK
CUTOFF 1000.0
UNITS k
NCONFIGS 2
ATOMS 1
LJ core 1.0  0.0
MOLTYPES 1
lj
MAXATOM 216
FINISH
VDW 1
LJ core  LJ core lj 1.0 1.0
CLOSE

The FIELD file one would use in a LSMC simulation is almost the same as one would use in any other simulation. The exception is that NCONFIGS must be set to 2 (see line 2 above). This is because two configurations are always used by DL_MONTE in a LSMC simulation. The reasons for this are technical, and explained in the user manual.

There are two somewhat esoteric aspects to the FIELD file above worth remarking upon. Firstly, we have set CUTOFF to 1000.0. This amounts to applying no cut-off to the Lennard-Jones potential - as we said we would do at the beginning of this tutorial. Secondly, we are using an energy unit of \(k_B\) K (where K is the kelvin). With this in mind, the Lennard-Jones parameter \(\epsilon\) (the first ‘1.0’ on line 12) is set to 1 \(k_B\) K. Moreover the parameter \(\sigma\) (the second ‘1.0’ on line 12) is set to 1 angstrom. This choice of units is for convenience. Recall that we wish to consider a temperature \(T=0.1\epsilon/k_B\). This temperature, given that \(\epsilon=1k_B\) K, corresponds to \(T=0.1\) K - which is reflected in the CONTROL file provided in a moment.

CONFIG, CONFIG.1 and CONFIG.2

We now consider the configuration files CONFIG, CONFIG.1 and CONFIG.2. As alluded to above, CONFIG is expected to contain two configurations in LSMC. These are prospective starting configurations for each phase. If the simulation is to be started in phase 1 (via a directive initactive in the CONTROL file - described in a moment), then the first configuration in CONFIG will be used as the starting configuration for the simulation. On the other hand if the simulation is to be started in phase 2, then the second configuration in CONFIG will be used as the starting configuration. By contrast CONFIG.1 and CONFIG.2 each contain only one configuration. The purpose of CONFIG.1 and CONFIG.2 is to define the underlying lattices to be used in lattice switch moves: the cell and particle positions in CONFIG.1 define the lattice for phase 1; and the cell and particle positions in CONFIG.2 define the lattice for phase 2. Details of exactly how the lattice switch move is constructed from the information in CONFIG.1 and CONFIG.2 can be found in the user manual.

Here we choose phase 1 to be hcp and phase 2 to be fcc. Thus CONFIG.1 and CONFIG.2 correspond to perfect (i.e. undistorted) hcp and fcc crystals respectively. Moreover we choose to start the simulation with the particle positions forming a perfect crystal. Thus the first configuration in CONFIG corresponds to a perfect hcp crystal, and the second configuration in CONFIG corresponds to a perfect fcc crystal. Of course one is free to use ‘equilibrated’ configurations, where the particle positions form an approximate crystal lattice, as prospective starting configurations in CONFIG.

CONTROL

Below is the CONTROL file. Note that in DL_MONTE lattice-switch Monte Carlo is referred to, for reasons we won’t discuss here, as phase-switch Monte Carlo (PSMC). You should regard the terms PSMC and LSMC as synonymous here.

PSMC simulation LJ                                      # Comment line

  use fed psmc                                          # Use PSMC in 'FED' functionality, followed by block of PSMC flags; begins FED block

      switchfreq 0                                      # Lattice-switch move frequency; NO SWITCH MOVES
      initactive 1                                      # Phase to start in (1 or 2; chooses configuration in CONFIG)

    psmc done                                           # End of PSMC block

    fed method ee 1.0 1.0 1000000000                    # Use EE method for learning the bias function; NEVER UPDATE DURING SIM; NO BIASING!
    fed order psmc 100 -1.0E10 1.0E10 1                 # Use PSMC order parameter for biasing; USE INFINITE RANGE

  fed done                                              # End of FED block

finish                                                  # End of 'use' block


ranseed                                                 # Use a random seed

steps 320000                                            # Simulation length

temperature 0.1
pressure 0.0

sample coordinate 21600                                 # Frequency to store configurations during simulation
archiveformat dlpoly4                                   # Configurations during simulation stored in DL_POLY-4 format (HISTORY)

yamldata 216                                            # Frequency to store energies etc. in yaml format (YAMLDATA)

move atoms 1 216                                        # Move atoms, followed by frequency
LJ core

move volume ortho log 1                                 # Move volume, followed by freqency


start

The important features of this file are as follows. Firsty, we have used the use fed psmc directive (line 3), which brings the LSMC machinery within DL_MONTE into play. This directive doubles as the start of the ‘PSMC block’ for specifying directives specific to LSMC, the end of which is signified by the psmc done directive (line 8). Here we see that there are two directives in this block, switchfreq and initactive. These are the only compulsory directives within the PSMC block; there are other optional directives for controlling various aspects of LSMC within DL_MONTE, but these are not required here. switchfreq specifies the frequency with which lattice switch moves are to be attempted. We have set this to 0 here because lattice switch moves are not required yet: recall that the goal of our preliminary simulations is to determine the range of order parameter exhibited by the equilibrium configurations of each phase, for which lattice switch moves are unnecessary (and, in fact, counterproductive: enabling lattice switch moves at this point would result in the system immediately jumping into the phase of lowest energy and remaining there during the preliminary simulation, as opposed to exploring the particular phase we are interested in). The other compulsory directive is initactive, which specifies which phase to start in; the corresponding configuration in CONFIG is used as the starting configuration.

Moving beyond the PSMC block, we come to the two compulsory FED directives, fed method and fed order [param] (lines 10 and 11). Note that in the latter we have specified, via fed order psmc, that the order parameter be the ‘PSMC’ order parameter. Moreover we have specified that the order parameter range to consider be split into 100 bins/states (via the first argument); but curiously we have chosen to consider an order parameter range between -1.0E10 and 1.0E10 (second and third arguments). In FED calculations moves which take the system outwith the specified range are rejected; the system is constrained to reside within the specified range. Since the aim of our preliminary simulations is to probe what would be a reasonable range to specify for later FED calculations, we do not wish to artificially constrain the range here, which is why we have set the range to be ‘infinite’. In a similar vein, while we choose in line 10 to use the expanded ensemble method (described later) to learn the bias function, we specify an update time larger than the length of the simulation (namely, 1,000,000,000 moves) to essentially ‘switch off’ this functionality for now, ensuring that the simulation is unbiased and thus visits the equilibrium configurations as would normally be the case. (Of course one does not necessarily have to use the expanded method at this point. Any method for learning the bias function would be fine to specify with fed method here - so long as biasing is essentially switched off).

The above ‘hacks’ may all seem a little strange, but it is currently the only way in DL_MONTE to track the LSMC order parameter while performing a ‘conventional’ Monte Carlo simulation, i.e. a Monte Carlo simulation without using lattice switch moves or biasing - as is required in our preliminary simulations.

Exercise

The input files correspond to a preliminary simulation for phase 1. Create a directory in which to perform this simulation, and copy the input files into it. Run the simulation. The simulation should take about 90 seconds to complete.

Now create another directory for the phase-2 preliminary simulation, and copy the input files into it. Change initactive in CONTROL so that the simulation is instead performed in phase 2, and run the simulation. Again, the simulation should take about 90 seconds to complete.

Post-processing YAMLDATA

Each preliminary simulation will output a file YAMLDATA.000, which contains a header of metadata (e.g. the temperature and pressure of the system), followed by data output periodically during the simulation (in YAML format). (The frequency of output to YAMLDATA.000 is controlled by the yamldata directive in CONTROL). Here are the first few lines of a YAMLDATA.000 output by an instance of the above phase-1 preliminary simulation:

temperature:       1.00000000000000E-01
pressure:          0.00000000000000E+00
usingbias:
usingpsmc:
---
-
  timestamp:                      216
  bias:          0.00000000000000E+00
  orderparam:    2.30319435732781E-01
  psmcactive:                       1
  energy:       -1.80527415326411E+03
  enthalpy:     -1.80527415326411E+03
  energyvdw:    -1.80527415326411E+03
  volume:        1.96415444935759E+02
-
  timestamp:                      432
  bias:          0.00000000000000E+00
  orderparam:    2.29806797369747E-01
  psmcactive:                       1
  energy:       -1.80506844807347E+03
  enthalpy:     -1.80506844807347E+03
  energyvdw:    -1.80506844807347E+03
  volume:        1.96357624881152E+02
-
  timestamp:                      648
  bias:          0.00000000000000E+00
  orderparam:    2.33436638488456E-01
  psmcactive:                       1
  energy:       -1.80424889431791E+03
  enthalpy:     -1.80424889431791E+03
  energyvdw:    -1.80424889431791E+03
  volume:        1.96387869892890E+02
-

The meaning of the data in the file should be transparent, though note that psmcactive: corresponds to the phase which the current configuration belongs to. We are interested in the range of \(M\) exhibited by the system in each phase, which corresponds to orderparam: in YAMLDATA.000. The script strip_yaml.sh in the script directory can be used to post-procedd YAMLDATA.000, extracting data from it and plotting it against simulation time (i.e. the number of moves so far). The argument to the script determines which piece of data to extract, as well as the name of the file to store the extracted data in. For example strip_yaml.sh energy would create a file energy.dat which contains the energy of the system vs. time.

Exercise

Use strip_yaml.sh to create files energy.dat and volume.dat for the preliminary simulations for each phase, which contain, respectively, the energy vs. simulation time and the volume vs. simulation time for the preliminary simulations. Plot the files, and check that the preliminary simulations are well-equilibrated within the simulation time.

Example plots of energy.dat and volume.dat from both preliminary simulations are given below.

Energy vs. simulation time for preliminary simulations. Note that the energy equilibrates quickly, within the first 50,000 moves. Note also that the energies of the hcp and fcc phases are very similar.

Volume vs. simulation time for preliminary simulations. Note that the volume equilibrates quickly, within the first 50,000 moves. Note also that the volumes of the hcp and fcc phases are very similar.

Use strip_yaml.sh to create files orderparam.dat for the preliminary simulations for each phase, which contain \(M\) vs. simulation time for the preliminary simulations. Plot the files and come up with a range of \(M\) which encompasses the equilibrium configurations for both phases.

An example plot of orderparam.dat for both preliminary simulations is given below.

Order parameter vs. simulation time for preliminary simulations. Note that the order parameter equilibrates quickly, within the first 50,000 moves. Note also that the order parameter distinguishes the hcp and fcc phases: the order parameter is negative for hcp configurations, while it is positive for fcc configurations.

Discussion

From the above figure it is clear that the hcp phase corresponds to an order parameter range of about -11 to -3, and the fcc phase corresponds to a range of about 2 to 11. Thus the range -11 to 11 encompasses both phases. However we will use the larger range -15 to 15 be safe. Configurations with order parameters less than -11 or greater than 11 may occur rarely at equilibrium though not be visible on our plot. However they still may make a significant contribution to the hcp-fcc free energy difference. However we expect that configurations with order parameters outwith the range -15 to 15 will be so rare to be insignificant.

Bias-generation simulation

Having determined an appropriate range for the order parameter (-15 to 15), we are now in a position to perform the bias-generation simulation.

CONTROL

For this simulation the files FIELD, CONFIG, CONFIG.1 and CONFIG.2 are identical to the preliminary simulations. However the CONTROL file is slightly different:

PSMC simulation LJ                                      # Comment line

  use fed psmc                                          # Use PSMC in 'FED' functionality, followed by block of PSMC flags; begins FED block

      switchfreq 216                                    # Lattice-switch move frequency
      initactive 1                                      # Phase to start in (1 or 2; chooses configuration in CONFIG)

      datafreq 43200                                    # Output frequency to PSDATA* (default = 1000)

    psmc done                                           # End of PSMC block

    fed method ee 1.0 1.0 43200000                      # Use EE method for learning the bias function; last argument is update frequency
    fed order psmc 100 -15.0 15.0 1                     # Bias over order parameter for PSMC in specified range

  fed done                                              # End of FED block

finish                                                  # End of 'use' block


ranseed                                                 # Use a random seed

steps 432000000                                         # Simulation length

temperature 0.1
pressure 0.0

sample coordinate 4320000                               # Frequency to store configurations during simulation
archiveformat dlpoly4                                   # Configurations during simulation stored in DL_POLY-4 format (HISTORY)

yamldata 43200                                          # Frequency to store energies etc. in yaml format (YAMLDATA)

move atoms 1 216                                        # Move atoms, followed by frequency
LJ core

move volume ortho log 1                                 # Move volume, followed by freqency

print 43200000                                          # Output frequency to OUTPUT.000 (not needed here - set frequency large)
stat  43200000                                          # Output frequency to PTFILE.000 (not needed here - set frequency large)

start

The important features of this file are as follows. Firstly, switchfreq (line 5) is no longer 0; lattice switch moves are enabled. Note that the frequency of lattice switch moves is high: here we attempt lattice switch moves with the same frequency as atom translation moves (i.e. with frequency 216). Because of the way LSMC is implemented in DL_MONTE, there is essentially no additional computational cost associated with attempting lattice switch moves. (Without going into detail, the computational cost of calculating the energy of the trial configuration generated by a lattice switch move is, in fact, absorbed into the calculation of the LSMC order parameter). Hence one may as well attempt lattice switch moves frequently, which ultimately improves the rate of successful transitions between phases.

Secondly, we have set the order parameter range to be from -15 to 15 (line 13). Moreover we have divided the order parameter range into 100 ‘states’ - the grid on which the bias function is defined.

Thirdly, we have now ‘enabled’ biasing (line 12), via the expanded ensemble method. This is the simplest method for learning the bias function. In the expanded ensemble method the bias is updated after every block of \(n\) moves. During each block we track how many the times the system visits each order parameter state. This information is then used to generate the bias function for the next block using the following equation:

\[\eta_i^{(k+1)}= \eta_i^{(k)}+\ln(h^{(k)}_i)\]

where \(\eta_i^{(k)}\) is the bias for state \(i\) during block \(k\), and \(h^{(k)}_i\) is the number of visits to state \(i\) during block \(k\). This equation increases the bias for states which were visited rarely during the last block, ensuring that they are visited more often during the next iteration. The end result is that after enough blocks all states are sampled equally. The forthcoming exercise will illustrate this. Here we have chosen \(n\) to be 43,200,000, which corresponds to every 100,000 ‘sweeps’ of the system (where sweep here is comprised with, on average, 216 atom translation moves, 216 lattice switch moves, and 1 volume move). This is the origin of the third argument to fed method ee on line 12. Note that it is necessary for \(n\) to be large in the expanded ensemble method - the block must be large enough that a representative histogram \(h^{(k)}_i\) is obtained for the current bias function \(\eta_i^{(k)}\). (Do not worry about the first two arguments to fed method ee, which are both ‘1.0’. By modifying these you can alter the form of the above equation. However this is unimportant here).

Fourthly, we have set the length of the simulation to be 432,000,000 moves; it can take a long time for an accurate bias function to be obtained. Note that this corresponds to 10 updates of the bias function.

Finally, though this is a minor point, we have used the sample coordinate, print and stat directives so that we only occasionally output data to the files ARCHIVE.000, OUTPUT.000 and PTFILE.000 - in order to keep the size of these files, which we do not use here, small. The directive datafreq in the PSMC block has the same purpose for the LSMC-specific output file PSDATA.000, which we also do not use here.

Exercise

Create a new directory in which to perform the bias-generation simulation. Copy the FIELD, CONFIG, CONFIG.1 and CONFIG.2 input files for one of the preliminary simulations into this directory, and create a new CONTROL file with the content given above. Run the simulation, noting that it should take about 5-6 hours to complete.

The first and third columns in the files FEDDAT.000_001, FEDDAT.000_002, …, FEDDAT.000_010 contain histograms over \(M\) (i.e. how often each order parameter state was visited) for blocks 1, 2, … 10 during the simulation. Plot the first and third columns in the FEDDAT file (e.g. using the command

plot "FEDDAT.000_001" using 1:3,"FEDDAT.000_002" using 1:3",...

in gnuplot) to inspect how the range of \(M\) explored during each block changes as the block number (i.e. the index on the FEDDAT file name) increases. You should see that the explored range increases with block number, until eventually the whole range is explored during a single block, as in the example plot below:

Order-parameter histograms for various blocks during a bias-generation simulation using the expanded ensemble method. Note that the histograms become increasingly flat as the block number is increased, and for blocks late in the simulation the entire order-parameter range is explored.

As you did previously during the preliminary simulations, use the strip_yaml.sh script to create a file orderparam.dat containing the order parameter vs. simulation time for the bias-generation simulation, and plot this. You should see that in the early stages of the simulation the order parameter is confined to small ranges of \(M\), but eventually the simulation sweeps quickly across the entire range of \(M\) - in a manner which reflects your previous plot. The orderparam.dat corresponding to the above plot is given below:

Order parameter vs. simulation time for a bias-generation simulation using the expanded ensemble method.

The first and second columns in the files FEDDAT.000_001, FEDDAT.000_002, …, FEDDAT.000_010 contain the bias functions calculated at the end of blocks 1, 2, …, 10. Plot the first two columns in the FEDDAT files to inspect how the bias function evolves throughout the simulation. You should see that it converges towards a double-well shape, as in the example plot below:

Evolution of the bias function during a bias-generation simulation using the expanded ensemble method.

Discussion

A well-converged bias function tells us a lot about the stability of the two phases in question; we need only actually proceed to the production simulation if we want quantitative results. It turns out that if the bias function is ‘perfect’, i.e. it leads to perfectly uniform sampling over the order parameter range, then the bias function is equivalent to the free energy profile over that order parameter. With this in mind, consider the final bias function obtained from the above simulation. Clearly it has two minima, one at \(M\approx -5\) and one at \(M\approx 5\). These are free energy minima corresponding to the equilibrium configurations of each phase: the \(M\approx -5\) corresponds to the hcp equilibrium configurations while the \(M\approx 5\) corresponds to the fcc equilibrium configurations. Between these minima is a maxima. This is the free energy barrier separating the two phases in the LSMC calcualtion. Note that the fcc minimum is lower than the hcp minimum, suggesting that the fcc phase is prefered at this temperature and pressure.

However, the above analysis supposes that the bias function obtained from the bias-generation simulation is perfect. It is not. Above we noted that the bias function is well converged. However it is not perfectly converged. Nevertheless at this point we can be fairly confident that the fcc phase is stable over hcp at this temperature and pressure. There is still the question though of what the precise value of the free energy difference is. We address this with the production simulation.

Bias-generation simulation (transition matrix)

This section is somewhat of an aside. Feel free to skip this section and proceed to the production simulation, perhaps returning to this section afterwards if you wish.

In ‘real world’ LSMC calculations the bias-generation simulation is the most technically difficult part. The bias function must be ‘good enough’ that it leads to both phases being explored in a reasonable timescale. However it becomes increasingly difficult to generate such a bias function as the system size increases. This is due to the fact that the height of the free energy barrier separating the two phases at \(M\approx 0\) increases with system size, making it increasingly difficult for the system to explore the entire range of \(M\) - even with biasing. It turns out that the method we used above, the expanded ensemble method, is only tractable for small systems. In this section we repeat the bias-generation simulation using a superior method, the transition matrix method, and contrast its performance against the expanded ensemble method.

The transition matrix method learns the bias function through tracking the proposed transitions between different order parameter states and their associated probabilities of being accepted. Details of the method can be found elsewhere (see the links at the end of this tutorial). As well as being more efficient (at least it is in this simulation) than the expanded ensemble method, the transition matrix method has the further benefit that it can can exploit parallelisation to reduce the computational time required to generate a good bias function. Multiple simulations can be run in parallel using the transition matrix method, and their results can be pooled together at their completion to generate a bias function. This is not possible with the expanded ensemble method. The simulations could even be assigned distinct regions of the considered order parameter range, reducing computational time further.

CONTROL

The CONTROL file for our bias-generation simulation exploiting the transition matrix method is as follows (all other input files are unchanged compared to the previous simulations):

PSMC simulation LJ                                      # Comment line

  use fed psmc                                          # Use PSMC in 'FED' functionality, followed by block of PSMC flags; begins FED block

      switchfreq 216                                    # Lattice-switch move frequency
      initactive 1                                      # Phase to start in (1 or 2; chooses configuration in CONFIG)

      datafreq 43200                                    # Output frequency to PSDATA* (default = 1000)

    psmc done                                           # End of PSMC block

    fed method tm 4320000 4320                          # Use TM method for learning the bias function; output and update frequencies
    fed order psmc 100 -15.0 15.0 1                     # Bias over order parameter for PSMC in specified range

  fed done                                              # End of FED block

finish                                                  # End of 'use' block

ranseed                                                 # Use a random seed

steps 43200000                                          # Simulation length

temperature 0.1
pressure 0.0

sample coordinate 4320000                               # Frequency to store configurations during simulation
archiveformat dlpoly4                                   # Configurations during simulation stored in DL_POLY-4 format (HISTORY)

yamldata 4320                                           # Frequency to store energies etc. in yaml format (YAMLDATA)

move atoms 1 216                                        # Move atoms, followed by frequency
LJ core

move volume ortho log 1                                 # Move volume, followed by freqency


print 43200000                                          # Output frequency to OUTPUT.000 (not needed here - set frequency large)
stat  43200000                                          # Output frequency to PTFILE.000 (not needed here - set frequency large)

start

Note that we have enabled the transition matrix method through the fed method tm directive on line 12. The first argument to this directive is how often the bias function is output to a FEDDAT file. Here we choose this to be every 4,320,000 moves. The second argument is how often the bias function is updated using the information regarding the transitions between order parameter states gathered so far in the simulation. It is beneficial to perform such updates frequenty, but it is unnecessary (and would result in an unnecessary computational expense) to do so every move. Here we set the updates to occur every 4320 moves, which corresponds to every 10 ‘sweeps’ of the system (where a sweep here is comprised with, on average, 216 atom translation moves, 216 lattice switch moves, and 1 volume move).

Note also that this simulation is far shorter than the previous bias-generation simulation: this simulation is of length 43,200,000 moves, while the previous one using the expanded ensemble method was of length 432,000,000 moves - a factor of 10 longer. We use such a short simulation here in anticipation of the fact that the transition matrix method is far more efficient than the expanded ensemble method, and hence will yield a good bias function in far less moves. As we shall see in a moment, 43,200,000 moves is indeed sufficient for the transition method to yield a good bias function, while this was far from the case using the expanded ensemble.

Exercise

Create a new directory in which to perform the bias-generation simulation using the transition matrix method. Copy the FIELD, CONFIG, CONFIG.1 and CONFIG.2 input files for one of the preliminary simulations into this directory, and create a new CONTROL file with the content given above. Run the simulation, noting that it should take about 30 minutes to complete.

As you did previously for the bias-generation simulation using the expanded ensemble method, plot the FEDDAT files and inspect the bias function throughout this bias-generation simulation. You should see that the bias function is well converged on the timescale of the simulation, as in the example plot below:

Evolution of the bias function during a bias-generation simulation using the transition matrix method.

Discussion

Note that here the transition matrix method converged to the correct bias function an order of magnitude faster than the expanded ensemble method. This serves to illustrate that significant speed-ups can be realised by using a more efficient algorithm. In fact, the transition matrix method can be exploited to realise even greater speed-ups than demonstrated here, by exploiting parallelisation. The relevant functionality exists in DL_MONTE. However it is beyond the scope of the current tutorial to demonstrate it.

Production simulation

We now turn to the production simulation, using the bias function obtained from the bias-generation simulation. Recall the aim of this simulation: to, with the aid of lattice switch moves and biasing, sample from both phases, and measure \(p_1/p_2\) (where \(p_1\) and \(p_2\) are the probabilities of the system being in phases 1 and 2 respectively), which can then be used to obtain the free energy difference \(\Delta G\) via the equation

\[\Delta G \equiv G_1 - G_2 = -k_BT\log(p_1/p_2)\]

There is, of course, the question of how we obtain \(p_1/p_2\) from the production simulation. After all, the production simulation uses biasing: how often phases 1 and 2 are visited during the production simulation does not reflect how often phases 1 and 2 are visited in the ‘true’ system (the quantities \(p_1\) and math:p_2 pertain to the true system). This is because we ‘force’ the simulation to visit both phases approximately equally using biasing. How often phases 1 and 2 are visited in the simulation in fact reflects the probabilities of the system being in either phase in the biased system. What we must do to obtain \(p_1\) and \(p_2\) (and hence \(p_1/p_2\)) is to take the data from the production simulation, and reweight it so that it pertains to the unbiased system. We will address this later, but first we must run the production simulation…

Input files

The input files for the production simulation are identical as for the bias-generation, with the following exceptions.

Firstly, we specify the bias function from the outset via an input file FEDDAT.000 - as opposed to learning the bias function during the simulation. The last FEDDAT file output from the bias-generation simulation, FEDDAT.000_010, contains the bias function we wish to use in the production simulation. Conveniently, this file can be used, without modification, as the input file FEDDAT.000 for the production simulation in order to specify the bias function. However one must instruct DL_MONTE that the bias function is to be read from FEDDAT.000, as opposed to being initialised as zero for all order parameter states as was the case in previous simulations. This is done by changing

fed order psmc 100 -15.0 15.0 1

in CONTROL to

fed order psmc 100 -15.0 15.0 -1

where the ‘-1’ signifies that the bias function is to be read from FEDDAT.000.

The other change which we must make to the input files, which is also to CONTROL, is that we must ensure that the bias function is not updated during the simulation, since we wish the bias function to remain fixed throughout. Recall that this was accomplished for the preliminary simulations by using the directive

fed method ee 1.0 1.0 1000000000

This directive is also suitable here.

Exercise

Create a new directory in which to perform the production simulation. Copy the FIELD, CONFIG, CONFIG.1, CONFIG.2 and CONTROL input files for the bias-generation simulation into this directory. Copy the FEDDAT.000_010 output file from the bias-generation simulation to this directory, and rename it FEDDAT.000. Modify the CONTROL file as described above so that it corresponds to a production simulation: this entails modifying the two lines containing the fed order and fed method directives as described above. Run the simulation, noting that it should take about 5-6 hours to complete.

Once the simulation is complete, there will be a file FEDDAT.000_001 which contains the bias function (first and second columns) and a histogram over \(M\) for the whole simulation (first and third columns). Note that the bias function will be identical to that in FEDDAT.000, since we instructed DL_MONTE not to update it during the simulation. As you did earlier for the FEDDAT files output by the bias-generation simulation, plot the first and third columns in FEDDAT.000_001 to see the aforementioned histogram. You should see that the entire order parameter range was sampled approximately uniformly, as was desired - as in the example plot below:

Order-parameter histogram for the production simulation.

As you did previously during the preliminary and bias-generation simulations, use the strip_yaml.sh script to create a file orderparam.dat containing the order parameter vs. simulation time, and plot this. You should see that the simulation efficiently explores the entire order parameter range, as in the example plot below:

Order parameter vs. simulation time for a production simulation.

Reweighting to remove the bias

Above you plotted the order parameter histogram for the production simulation. This information was contained in the first and third columns of the file FEDDAT.000_001. Let \(h_i\) denote how often order parameter state \(i\) was visited during the simulation. Conveniently, the bias \(\eta_i\) for each order parameter state \(i\) is contained in the second column of FEDDAT.000_001. Now, the effect of biasing is to add an energy \(-\eta_i/\beta\) to all configurations belonging to order parameter state \(i\). This amplifies the sampling of state \(i\) by a factor \(\exp(\eta_i)\): if \(i\) were visited with frequency \(f\) during an unbiased simulation, it would be visited with frequency \(f\exp(\eta_i)\) in a biased simulation. Conversely if \(i\) were visited with frequency \(f'\) in a biased simulation, it would be visited with frequency \(f'\exp(-\eta_i)\) in an (ergodic) unbiased simulation. Thus we can obtain the order-parameter histogram corresponding to an unbiased simulation by multiplying the frequency corresponding to each state in the biased simulation by \(\exp(-\eta_i)\).

Exercise

The script unbiasedhist.sh performs the aforementioned reweighting, using the data in FEDDAT.000_001 file in the current directory. It creates a file unbiasedhist.dat which contains the corresponding unbiased histogram. Apply this script to the output of your production simulation, and plot the resulting file unbiasedhist.dat. An example plot is given below:

Unbiased order-parameter histogram created by reweighting the biased order-parameter histogram obtained from the production simulation.

Note that there are two peaks in the unbiased histogram. These correspond to the equilibrium configurations for hcp (left) and fcc (right) phases. Note also that the fcc peak is larger than the hcp peak, signifying that the fcc phase is stable over hcp at this temperature and pressure. Thus what we saw previously earlier for the bias function while analysing our data from the bias-generation simulation is borne out in the data from the production simulation.

Recall that negative order parameters correspond to phase 1 while positive order parameters correspond to phase 2. With this in mind, the unbiased histogram can be used to obtain \(p_1/p_2\) as follows: sum all elements of the histogram corersponding to phase 1, do the same for elements corresponding to phase 2, and then divide the former by the latter. The script unbiasedhist_probratio.sh performs this task, acting on the unbiasedhist.dat file in the current directory. It outputs a single value, \(p_1/p_2\).

Apply unbiasedhist_probratio.sh to extract \(p_1/p_2\) for your production simulation. Substitute this quantity into the expression for \(\Delta G\) given at the beginning of this section in order to obtain a value for \(\Delta G\), remembering that our simulations used \(T=0.1\epsilon/k_B\). Divide this value of \(\Delta G\) by the number of particles in the system, \(N\), to get the free energy difference per particle \(\Delta G/N\) (in units of \(\epsilon\)).

The ‘true’ value of \(\Delta G/N\) for this system, obtained using long DL_MONTE simulations and verified against another LSMC code, is \(\Delta G/N=0.001283(7)\epsilon\), where the ‘7’ in parenthesis is the uncertainty on the final digit (the standard error of the mean - see below). Is your value of \(\Delta G/N\) in agreement with the true value?

Calculating uncertainties

Above we evaluated \(\Delta G\) by deducing an unbiased order-parameter histogram for the simulation, integrating over the regions of the histogram corresponding to phases 1 and 2 to obtain \(p_1/p_2\), and then substituting this \(p_1/p_2\) into the equation relating \(p_1/p_2\) to \(\Delta G\). This procedure certainly yields a value for \(\Delta G\). However there is the question of what the uncertainty in the \(\Delta G\) obtained from our production simulation is. The above method for obtaining \(\Delta G\) does not provide this. Moreover the histogram used in the analysis above takes into account the configurations visited at the very start of the simulation - during the equilibration time. These configurations should not be considered when evaluating \(\Delta G\) (or indeed any other quantity). Here we demonstrate how to calculate \(\Delta G\) and its associated uncertainty more rigourously using the data output by the production simulation.

We first briefly revise how to combine multiple measurements of a physical quantity \(X\) into a single ‘quoted value’ and an associated uncertainty. Let us assume that we have measured \(X\) \(n\) times. Let \(X_1,X_2,\dotsc,X_n\) denote these \(n\) measurements. The standard practice is to calculate the mean of all the measurements,

\[\bar{X}=\frac{1}{n}\sum_{i=1}^nX_i,\]

and then use this as the ‘quoted value’ of \(X\). Regarding the uncertainty in \(\bar{X}\), one commonly uses the standard error of the mean:

\[\delta X=\text{stdev}(X)/\sqrt{n}\]

as the uncertainty in \(\bar{X}\), where \(\text{stdev}(X)\) is the standard deviation of the observations \(X_1,X_2,\dotsc,X_n\). However it is also common to use twice the standard error of the mean as the uncertainty. (This corresponds to a 95% confidence interval, while the standard error of the mean corresponds to a 68% confidence interval).

Due to the use of random numbers in Monte Carlo, diffferent simulations (using different random-number-generator seeds) will yield different values for a physical quantity \(X\). We can treat the results of different simulations as different ‘measurements’ in the sense above, and combine the results from these simulations to obtain a quoted value and an associated uncertainty as described above. In fact, one does not need to resort to combining the results from separate simulations. Instead, one can split the data obtained from a single simulation into blocks, calculate \(X\) for each of these blocks, and use these \(X\) as the measurements from which the quoted value and uncertainty are calculated. We apply this approach below to the data output by the production simulation to calculate a quoted value and associated uncertainty for \(\Delta G\).

Exercise

The script probratio.sh processes the data in YAMLDATA.000 to calculate \(p_1/p_2\), but only the data in YAMLDATA.000 which corresponds to a specified ‘time window’. The script takes two arguments, \(t_{\text{min}}\) and \(t_{\text{max}}\), which correspond to the lower and upper bounds of the window - where the units of time here are the number of elapsed Monte Carlo moves (which corresponds to the timestamp: data in YAMLDATA.000). An example usage is as follows, where the script is called in the directory containing the YAMLDATA.000 file for the production simulation (the $ is the shell prompt):

$ probratio.sh 4320000 432000000
0.0610501

In this example \(t_{\text{min}}=4320000\) and \(t_{\text{max}}=432000000\), and 0.0610501 is the value of \(p_1/p_2\) corresponding to the this time window.

Recall that the length of our production simulation was 432,000,000 moves. This can be split into 4 blocks of 108,000,000 moves each, in which case probratio.sh can be used to obtain \(p_1/p_2\) obtained for each block via

probratio.sh 1 108000000

for the 1st block,

probratio.sh 108000001 216000000

for the 2nd block,

probratio.sh 216000001 324000000

for the 3rd block, and

probratio.sh 324000001 432000000

for the 4th block. The \(p_1/p_2\) obtained for each block in this manner can then be used to obtain a \(\Delta G\) for each block as described earlier, after which the mean and standard error of the mean of all blocks can be calculated. Do this for your production simulation, and compare your result to the ‘true’ value given above.

Do the same again, but this time assume an equilibration time of 4,320,000 moves, i.e. ignore the first 4,320,000 moves, partitioning the remaining time into four blocks.

One concern is how large each block should be. The answer is that the ‘block size’, i.e. the size of the time window each block coversm should be larger than the correlation time for the physical quantity \(X\) one is interested in. This ensures that the blocks, and hence also their associated values of \(X\), are uncorrelated. For the purpose of evaluating \(\Delta G\) the block size should be larger than the time it takes for the system to explore both phases. One can check this by plotting the order parameter against simulation time, and inspecting the timescale of the fluctuations in the order parameter against the block size: the timescale of the fluctuations should be far smaller than the block size. Verify that this is indeed the case here. (Recall that you have already plotted the order parameter against simulation time for the production simulation earlier).

As one final check, repeat the above for various block sizes. The essential result should be independent of block size - so long as one does not make the block size so small that both phases are not explored in any block, or that the block size is so small that time resolution of the data in YAMLDATA.000 is too coarse to yield meaningful results.

Further information

• Section 2 of the paper at https://arxiv.org/abs/1609.04329 (which describes an old LSMC code) provides a more thorough discussion of LSMC relevant to this tutorial, including the expanded ensemble method (called the ‘visited states method’ there), the transition matrix method, and a description of the LSMC order parameter \(M\) relevant to the simulations performed here
• The DL_MONTE reference manual provides more detailed information regarding LSMC functionality in DL_MONTE